##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)
library(ggpubr)
library(ggsci)
library(ggrepel)

##########################################################################################

option_list <- list(
    make_option(c("--input_file"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    input_file <- paste(work_dir,"/images/evolutionTime/CompareTrunkSMG.Time.tsv",sep="")
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.tsv"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/evolutionTime"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_file <- opt$input_file
sample_info <- opt$sample_info
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

dat_input <- data.frame(fread(input_file))
dat_info <- data.frame(fread(sample_info))

###########################################################################################

dat_info <- dat_info %>%
group_by( Patient , Age , Gender , Tobacco , Alcohol , PickleFood , HP , TCGA_Class ) %>%
summarize( BurdenAll = median(BurdenAll) ,BurdenExon = median(BurdenExon) )

dat_info <- unique(subset( dat_info , TCGA_Class != "IM" ))[,c(1:8)]

###########################################################################################

dat_input2 <- merge( dat_input , dat_info , by.x = "ID" , by.y = "Patient" )
dat_input2$Class <- factor( dat_input2$Class , levels = c("IGC" , "DGC") , order = T )

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("P == ",down," %*% 10","^",up)
    return(text)
}

col_tmp <- c(
    rgb(red=179,green=60,blue=59,alpha=255,max=255) ,
    rgb(red=14,green=90,blue=170,alpha=255,max=255)
)

y_max <- max(dat_input2$mean) + 1
y_lab <- "Number of years for \nIM to GC progression"

base_col <- c("Gender" , "Tobacco" , "Alcohol" , "HP" , "TCGA_Class")

dat_input2_all <- dat_input2
dat_input2_all$Class <- "All"

dat_input3 <- rbind( dat_input2 , dat_input2_all )
dat_input3$Class <- factor( dat_input3$Class , levels = c( "All" , "IGC" , "DGC") , order = T )
dat_input3$TCGA_Class <- factor( dat_input3$TCGA_Class , levels = c("GS" , "CIN") , order = T )

##################################################################################

for( baseuse in base_col ){

    dat_plot_tmp <- data.frame(dat_input3)
    dat_plot_tmp$useCol <- dat_plot_tmp[[baseuse]]
    dat_plot_tmp <- subset( dat_plot_tmp , !is.na(useCol) )

    dat_tmp <- c()
    for(class in unique(dat_plot_tmp$Class) ){

        dat_plot_tmp_use <- subset( dat_plot_tmp , Class == class )

        p <- wilcox.test(
            subset(dat_plot_tmp_use , useCol==unique(dat_plot_tmp_use$useCol)[1])$mean , 
            subset(dat_plot_tmp_use , useCol==unique(dat_plot_tmp_use$useCol)[2])$mean)$p.value

        if( p < 0.01 ){
            p_text <- trans(p)
        }else{
            p_text <- paste0( "P == " , round(as.numeric(p) , 3) ) 
        }

        dat_plot_tmp_use$p_text <- ""
        dat_plot_tmp_use$p_text[1] <- p_text

        dat_tmp <- rbind( dat_plot_tmp_use , dat_tmp )
    }


    plot <- ggplot( dat_tmp , aes( x = useCol , y = mean , color = useCol ) ) +
        geom_boxplot(size = 1.2 , outlier.shape = NA ) + ## 去除散点，加粗线
        scale_color_manual(values=col_tmp) +
        scale_fill_manual(values =col_tmp) +
        geom_jitter(position=position_jitter(0.2 , seed = 1),aes(color=useCol)) +
        geom_text(aes(label=p_text , y = y_max , x = 1.5),parse = TRUE,size=5 , color = "black", face='bold') +
        xlab(NULL) +
        ylab(y_lab)+
        theme_bw() +
        facet_grid(.~Class) +
        theme(
            legend.position = 'none',
            legend.title = element_blank() ,
            panel.grid.major=element_blank(),
            panel.grid.minor=element_blank(),
            panel.background = element_blank(),
            panel.border = element_blank(),
            plot.title = element_text(size = 12,color="black",face='bold'),
            legend.text = element_text(size = 12,color="black",face='bold'),
            axis.text.y = element_text(size = 12,color="black",face='bold'),
            axis.title.x = element_text(size = 12,color="black",face='bold'),
            axis.title.y = element_text(size = 14,color="black",face='bold'),
            #axis.text.x = element_text(size = 14,color="black",face='bold',angle = 45,hjust = 1) ,
            axis.text.x = element_text(size = 14,color="black",face='bold') ,
            axis.ticks.length = unit(0.2, "cm") ,
            strip.text.x = element_text(size = 17, colour = "black",face='bold') ,
            axis.line = element_line(size = 0.5)) 

    width <- 5/1
    height <- 4.47/1
    out_name <- paste0(out_path , "/CompareTrunkSMG.Time." , baseuse , ".pdf")
    ggsave(file=out_name,plot=plot,width=width,height=height)
}